
#include "rbfFunction_3to1.h"
#include "rbfFunction_3to3.h"
#include "rbfEvalKernel_3to1.h"
#include "rbfEvalKernel_3to3.h"

#include "rbfUtils.h"
#include "typedefs.h"
#include "MCpolygonizer.h"

#include <iostream>
#include <cmath>
#include <iterator>
#include <numeric>
#include <vector>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>

using namespace std;

const int NCENTERSLIM=16;
const int NPOINTSLIM=16;

/// Returns {max, min, avg, perOpAvg} error values given two arrays of data
void errorValues(const vector<Number> &, const vector<Number> &, Number []);
/// Returns {max, min, avg, perOpAvg} error values given two arrays of data
void errorValues(const vector<Point3> &, const vector<Point3> &, Number []);


void runRbfEval_3to1_randomData_CPU();
void runRbfEval_3to1_randomData_GPUfs();
void runRbfEval_3to1_randomData_error();

void runRbfEval_3to3_randomData_CPU();
void runRbfEval_3to3_randomData_GPUfs();
void runRbfEval_3to3_randomData_error();

void runRbfEval_3to1();
void runRbfEval_3to3();

int main(int argc, char **argv) {
   
   glutInit(&argc, argv);
   GLuint winHandle = glutCreateWindow(argv[0]);
   
   //~ runRbfEval_3to1_randomData_CPU();
   //~ runRbfEval_3to1_randomData_GPUfs();
   runRbfEval_3to1_randomData_error();
   
   //~ runRbfEval_3to3_randomData_CPU();
   //~ runRbfEval_3to3_randomData_GPUfs();
   //~ runRbfEval_3to3_randomData_error();
   
   //~ runRbfEval_3to1();
   //~ runRbfEval_3to3();
   
   glutDestroyWindow(winHandle);
   return 0;
}

void runRbfEval_3to1_randomData_CPU() {
   
   cout<<endl<<"Biharmonic splines (3to1) evaluation on CPU ..."<<endl;
   time_t startEx = clock();
   for (int i=12; i<=NCENTERSLIM; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Number> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), randomGenerator);
      
      RbfFunction3to1 rbfeval;
      rbfeval.setCentersAndWeights(centers, weights);
      for (int j=12; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         
         rbfeval.setEvalPoints(points);
         vector<Number> fvalues(points.size());
         start = clock();
         rbfeval.execute(fvalues);
         end = clock();
         Number totaltime = (end-start)/(double)CLOCKS_PER_SEC;
         cout<<setw(7)<<totaltime<<"  |  ";
         
         /// file output
         //~ ostringstream nss, mss; nss<<i; mss<<j;
         //~ string filename = "cpuRbfEval3to1_2(" + mss.str() + ")x_2(" + nss.str() + ").txt";
         //~ ofstream fileout(filename.c_str());
         //~ fileout.precision(8);
         //~ fileout<<"cx cy cz w"<<endl;
         //~ for (int ii=0; ii<centers.size(); ++ii) fileout<<centers[ii]<<" "<<weights[ii]<<endl;
         //~ fileout<<"px py pz fvalue"<<endl;
         //~ for (int ii=0; ii<points.size(); ++ii)  fileout<<points[ii]<<" "<<fvalues[ii]<<endl;
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

void runRbfEval_3to1_randomData_GPUfs() {
   
   cout<<endl<<"Biharmonic splines (3to1) evaluation on GPU (fragment shader)..."<<endl;
   time_t startEx = clock();
   for (int i=12; i<=NCENTERSLIM; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Number> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), randomGenerator);
      
      RbfEvalKernel3to1 rbfeval;
      rbfeval.setCentersAndWeights(centers, weights);
      for (int j=12; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         
         rbfeval.setEvalPoints(points);
         int N = (int)(ceil(sqrt(points.size())));
         vector<Number> fvalues(N*N); // set the correct size
         start = clock();
         rbfeval.execute(fvalues);
         end = clock();
         Number totaltime = (end-start)/(double)CLOCKS_PER_SEC;
         cout<<setw(7)<<totaltime<<"  |  ";
         
         /// file output
         //~ ostringstream nss, mss; nss<<i; mss<<j;
         //~ string filename = "gpuRbfEval3to1_2(" + mss.str() + ")x_2(" + nss.str() + ").txt";
         //~ ofstream fileout(filename.c_str());
         //~ fileout.precision(8);
         //~ fileout<<"cx cy cz w"<<endl;
         //~ for (int ii=0; ii<centers.size(); ++ii) fileout<<centers[ii]<<" "<<weights[ii]<<endl;
         //~ fileout<<"px py pz fvalue"<<endl;
         //~ for (int ii=0; ii<points.size(); ++ii)  fileout<<points[ii]<<" "<<fvalues[ii]<<endl;
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

void runRbfEval_3to1_randomData_error() {
   
   cout<<endl<<"Biharmonic splines (3to1) evaluation (error comparison)..."<<endl;
   time_t startEx = clock();
   //~ for (int i=16; i<=NCENTERSLIM; i+=2) {
   for (int i=15; i<=15; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Number> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), randomGenerator);
      
      RbfEvalKernel3to1 rbfeval_fs; // The RBF 
      rbfeval_fs.setCentersAndWeights(centers, weights);
      
      RbfFunction3to1 rbfeval_cpu;
      rbfeval_cpu.setCentersAndWeights(centers, weights);
      
      for (int j=16; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         
         // CPU
         rbfeval_cpu.setEvalPoints(points);
         vector<Number> fvalues_cpu(points.size());
         rbfeval_cpu.execute(fvalues_cpu);
         
         // GPU (fragment shader)
         rbfeval_fs.setEvalPoints(points);
         int N = (int)(ceil(sqrt(points.size())));
         vector<Number> fvalues_fs(N*N); // set the correct size
         rbfeval_fs.execute(fvalues_fs);
         
         Number error[4];
         errorValues(fvalues_cpu, fvalues_fs, error);
         
         // error min, error max, average error
         //~ cout<<setw(7)<<error[0]<<" : "<<error[1]<<" : "<<error[2]<<" : "<<error[3]<<"  ||  ";
         cout<<setw(7)<<error[0]<<"  ||  ";
         
         /// file output
         ostringstream nss, mss; nss<<i; mss<<j;
         string filename = "comparisonEval3to1_2(" + mss.str() + ")x_2(" + nss.str() + ").txt";
         ofstream fileout(filename.c_str());
         fileout.precision(8);
         fileout<<"px py pz fvalues_cpu fvalues_fs"<<endl;
         for (int ii=0; ii<points.size(); ++ii)  fileout<<points[ii]<<" "<<fvalues_cpu[ii]<<" :: "<<fvalues_fs[ii]<<endl;
         fileout<<"cx cy cz w"<<endl;
         for (int ii=0; ii<centers.size(); ++ii) fileout<<centers[ii]<<" "<<weights[ii]<<endl;
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

/// ////////////////////////////

void runRbfEval_3to3_randomData_CPU() {
   
   cout<<endl<<"Biharmonic splines (3to3) evaluation on CPU ..."<<endl;
   time_t startEx = clock();
   for (int i=12; i<=NCENTERSLIM; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Point3> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), pointRandomGenerator);
      
      RbfFunction3to3 rbfeval;
      rbfeval.setCentersAndWeights(centers, weights);
      for (int j=12; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         rbfeval.setEvalPoints(points);
         vector<Point3> fvalues(points.size());
         start = clock();
         rbfeval.execute(fvalues);
         end = clock();
         Number totaltime = (end-start)/(double)CLOCKS_PER_SEC;
         cout<<setw(7)<<totaltime<<"  |  ";
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

void runRbfEval_3to3_randomData_GPUfs() {
   
   cout<<endl<<"Biharmonic splines (3to3) evaluation on GPU (fragment shader) ..."<<endl;
   time_t startEx = clock();
   for (int i=12; i<=NCENTERSLIM; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Point3> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), pointRandomGenerator);
      
      RbfEvalKernel3to3 rbfeval;
      rbfeval.setCentersAndWeights(centers, weights);
      for (int j=12; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         rbfeval.setEvalPoints(points);
         
         int N = (int)(ceil(sqrt(points.size())));
         vector<Point3> fvalues(N*N); // set the correct size
         start = clock();
         rbfeval.execute(fvalues);
         end = clock();
         Number totaltime = (end-start)/(double)CLOCKS_PER_SEC;
         cout<<setw(7)<<totaltime<<"  |  ";
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

void runRbfEval_3to3_randomData_error() {
   
   cout<<endl<<"Biharmonic splines (3to3) evaluation (error comparison)..."<<endl;
   time_t startEx = clock();
   for (int i=12; i<=NCENTERSLIM; i+=2) {
      time_t start, end;
      int n = (int)pow(2.0,i);     // the number of centers
      vector<Point3> centers(n);    // the RBF's centers
      vector<Point3> weights(n);    // the RBF's weights
      generate(centers.begin(), centers.end(), pointRandomGenerator);
      generate(weights.begin(), weights.end(), pointRandomGenerator);
      
      RbfEvalKernel3to3 rbfeval_fs; // The RBF 
      rbfeval_fs.setCentersAndWeights(centers, weights);
      
      RbfFunction3to3 rbfeval_cpu;
      rbfeval_cpu.setCentersAndWeights(centers, weights);
      for (int j=12; j<=NPOINTSLIM; j+=2) {
         int m = (int)pow(2.0,j);  // the number of points
         vector<Point3> points(m);  // the points
         generate(points.begin(), points.end(), pointRandomGenerator);
         
         // CPU
         rbfeval_cpu.setEvalPoints(points);
         vector<Point3> fvalues_cpu(points.size());
         rbfeval_cpu.execute(fvalues_cpu);
         
         // GPU (fragment shader)
         rbfeval_fs.setEvalPoints(points);
         int N = (int)(ceil(sqrt(points.size())));
         vector<Point3> fvalues_fs(N*N); // set the correct size
         rbfeval_fs.execute(fvalues_fs);
         
         Number error[4];
         errorValues(fvalues_cpu, fvalues_fs, error);
         
         // error min, error max, average error
         //~ cout<<setw(7)<<error[0]<<" : "<<error[1]<<" : "<<error[2]<<" : "<<error[3]<<"  ||  ";
         cout<<setw(7)<<error[0]<<"  ||  ";
      }
      cout<<endl;
   }
   cout<<"Experiment total time: "<<(clock()-startEx)/(double)CLOCKS_PER_SEC<<"s."<<endl;
}

/// /////////////////////////////

void runRbfEval_3to1() {
   
   int nvGrid = 55; // number of vertices of grid
   time_t start, end;
   // the RBF's centers
   int ncenters;
   ifstream inputfile("knot-512v.vis");
   assert(inputfile);
   inputfile>>ncenters;
   vector<Point3> centers(ncenters);
   vector<Number> weights(ncenters);
   cout<<"centers number: "<<ncenters<<endl;
   for (int i=0; i<ncenters; ++i) {
      Point3 c; Number w;
      inputfile>>c[0]>>c[1]>>c[2]>>w;
      centers[i] = c;
      weights[i] = w;
   }
   assert(ncenters == centers.size());
   
   // the points to be evaluated
   vector<Point3> points;
   gridGenerator(nvGrid, nvGrid, nvGrid, points);
   
   /// /////////////////////
   cout<<"evaluation begin (CPU)..."<<endl;
   start = clock();
   RbfFunction3to1 rbf;
   rbf.setCentersAndWeights(centers, weights);
   rbf.setEvalPoints(points);
   vector<Number> fvalues(points.size());
   rbf.execute(fvalues);
   end = clock();
   std::cout<<"::: evaluation time (CPU): "<<(end-start)/(double)CLOCKS_PER_SEC<<"s."<<endl;
   cout<<"evaluation end"<<endl;
   
   {
      
   cout<<"evaluation begin (GPU)..."<<endl;
   start = clock();
   /// Kernel execution
   RbfEvalKernel3to1 kernel;
   kernel.setCentersAndWeights(centers, weights);
   kernel.setEvalPoints(points);
   
   int N = (int)(ceil(sqrt(points.size())));
   vector<Number> fvalues(N*N);
   kernel.execute(fvalues);
   
   //copy(result.begin(), result.end(), ostream_iterator<Center>(cout, " "));
   end = clock();
   //~ cout<<"points.size(): "<<points.size()<<endl;
   //~ cout<<"result.size(): "<<result.size()<<endl;
   //~ result.resize(points.size());
   //copy(result.begin(), result.end(), ostream_iterator<Center>(cout, " "));
   std::cout<<"::: evaluation total time GPU = "<<(end-start)/(double)CLOCKS_PER_SEC<<"s."<<endl;
   cout<<"evaluation end"<<endl;
   
   
   cout<<"polygonizer begin ..."<<endl;
   start = clock();
   MCPolygonizer poly(nvGrid, nvGrid, nvGrid);
   poly.polygonize(points, fvalues);
   poly.export2off("knot55.off");
   end = clock();
   std::cout<<"::isosurface extraction total time:: "<<(end-start)/(double)CLOCKS_PER_SEC<<"s."<<endl;
   cout<<"polygonizer end"<<endl;
   }
   
}

void runRbfEval_3to3() {
   
   time_t start, end;
   // the RBF's centers
   int ncenters = 6;
   vector<Point3> centers(ncenters);
   //~ generate(centers.begin(), centers.end(), pointRandomGenerator);
   centers[0] = Point3(0.5, -0.5, -0.5); // 1
   centers[1] = Point3(-0.5, -0.5, -0.5); // 2
   centers[2] = Point3(-0.5, -0.5, 0.5); // 6
   centers[3] = Point3(0.5, -0.5, 0.5); // 5
   centers[4] = Point3(0.5, 0.5, 0.5); // 4
   centers[5] = Point3(-0.5, 0.5, -0.5); // 3
   //~ cout<<"centers:"<<endl;
   //~ copy(centers.begin(), centers.end(), ostream_iterator<Point3>(cout, "\n"));
   
   // The points to be evaluated
   vector<Point3> points(8);
   points[0] = Point3(0.5, 0.5, -0.5); // 0
   points[1] = Point3(0.5, -0.5, -0.5); // 1
   points[2] = Point3(-0.5, -0.5, -0.5); // 2
   points[3] = Point3(-0.5, 0.5, -0.5); // 3
   points[4] = Point3(0.5, 0.5, 0.5); // 4
   points[5] = Point3(0.5, -0.5, 0.5); // 5
   points[6] = Point3(-0.5, -0.5, 0.5); // 6
   points[7] = Point3(-0.5, 0.5, 0.5); // 7
   //~ cout<<"points:"<<endl;
   //~ copy(points.begin(), points.end(), ostream_iterator<Point3>(cout, "\n"));
   //~ cout<<endl;
   
   cout<<"begin (CPU)..."<<endl;
   RbfFunction3to3 rbfcpu;
   rbfcpu.compute_iM(centers);
   //~ cout<<rbfcpu.iM<<endl<<endl;
   rbfcpu.setEvalPoints(points);
   for (int i=0; i<1; ++i) {
      vector<Point3> fcenters(ncenters);
      //~ generate(centersEnd.begin(), centersEnd.end(), pointRandomGenerator);
      fcenters[0] = Point3(0.5, -0.5, -0.5); // 1
      fcenters[1] = Point3(-0.5, -0.5, -0.5); // 2
      fcenters[2] = Point3(-0.5, -0.5, 0.5); // 6
      fcenters[3] = Point3(0.5, -0.5, 0.5); // 5
      fcenters[4] = Point3(0.2, 0.3, 0.45); // 4
      fcenters[5] = Point3(-0.5, 0.5, -0.5); // 3
      //~ cout<<"fcenters:"<<endl;
      //~ copy(fcenters.begin(), fcenters.end(), ostream_iterator<Point3>(cout, "\n"));
      rbfcpu.genWeights(fcenters);
      //~ cout<<"rbfcpu.weights:"<<endl;
      //~ copy(rbfcpu.weights.begin(), rbfcpu.weights.end(), ostream_iterator<Point3>(cout, "\n"));
      
      vector<Point3> result(points.size());
      start = clock();
      rbfcpu.execute(result);
      end = clock();
      result.resize(points.size());
      std::cout<<"::: total time CPU = "<<(end-start)/(double)CLOCKS_PER_SEC<<"s."<<endl;
      cout<<"result:"<<endl;
      copy(result.begin(), result.end(), ostream_iterator<Point3>(cout, "\n"));
   }
   cout<<"end (CPU)"<<endl;
   
   cout<<endl<<"::::::::::"<<endl<<endl;
   
   RbfEvalKernel3to3 kernel;
   kernel.compute_iM(centers);
   //~ cout<<kernel.iM<<endl<<endl;
   //~ copy(kernel.centers.begin(), kernel.centers.end(), ostream_iterator<Point3>(cout, "\n"));
   kernel.setEvalPoints(points);
   //~ copy(kernel.pPoints->begin(), kernel.pPoints->end(), ostream_iterator<Point3>(cout, "\n"));
   //~ kernel.install(true);
   for (int i=0; i<1; ++i) {
      // the RBF modified centers
      vector<Point3> fcenters(ncenters);
      //~ generate(centersEnd.begin(), centersEnd.end(), pointRandomGenerator);
      fcenters[0] = Point3(0.5, -0.5, -0.5); // 1
      fcenters[1] = Point3(-0.5, -0.5, -0.5); // 2
      fcenters[2] = Point3(-0.5, -0.5, 0.5); // 6
      fcenters[3] = Point3(0.5, -0.5, 0.5); // 5
      fcenters[4] = Point3(0.2, 0.3, 0.45); // 4
      fcenters[5] = Point3(-0.5, 0.5, -0.5); // 3
      
      kernel.genWeights(fcenters);
      //~ cout<<"weights GPU: "<<endl;
      //~ copy(kernel.weights.begin(), kernel.weights.end(), ostream_iterator<Point3>(cout, "\n"));
      
      cout<<"evaluation begin (GPU)..."<<endl;
      start = clock();
      
      /// Download the data
      int N = (int)(ceil(sqrt(points.size())));
      vector<Point3> result(N*N);
      /// Kernel execution
      kernel.execute(result);
      cout<<"result before"<<endl;
      //~ copy(result.begin(), result.end(), ostream_iterator<Point3>(cout, "\n"));
      kernel.addPolynomialTerm(result);
      end = clock();
      cout<<"result after"<<endl;
      result.resize(points.size());
      copy(result.begin(), result.end(), ostream_iterator<Point3>(cout, "\n"));
      std::cout<<"::: total time GPU = "<<(end-start)/(double)CLOCKS_PER_SEC<<"s."<<endl;
      cout<<"evaluation end"<<endl;
   }
   
}

/// Returns {max, min, avg, perOpAvg} error values given two arrays of data
void errorValues(const vector<Number> &fvalcpu, const vector<Number> &gpugrid, Number errorVal[]) {
   
   int npoints = fvalcpu.size();
   Number errorMax = 0, errorMin = RAND_MAX, errorAvg = 0, errorPerOpAvg = 0;
   for(int i=0; i<npoints; i++) {
      Number error = abs(fvalcpu[i] - gpugrid[i]);
      errorAvg += error;
      errorMax = max(error, errorMax);
      errorMin = min(error, errorMin);
   }
   
   errorAvg /= npoints;
   errorPerOpAvg = errorAvg/npoints;
   
   errorVal[0] = errorMax;
   errorVal[1] = errorMin;
   errorVal[2] = errorAvg;
   errorVal[3] = errorPerOpAvg;
}

/// Returns {max, min, avg, perOpAvg} error values given two arrays of data
void errorValues(const vector<Point3> &fvalcpu, const vector<Point3> &gpugrid, Number errorVal[]) {
   
   int npoints = fvalcpu.size();
   Number errorMax = 0, errorMin = RAND_MAX, errorAvg = 0, errorPerOpAvg = 0;
   for(int i=0; i<npoints; i++) {
      Vector3 v = fvalcpu[i] - gpugrid[i];
      Number error = sqrt(v*v); // the lenght of the difference vector
      errorAvg += error;
      errorMax = max(error, errorMax);
      errorMin = min(error, errorMin);
   }
   
   errorAvg /= npoints;
   errorPerOpAvg = errorAvg/npoints;
   
   errorVal[0] = errorMax;
   errorVal[1] = errorMin;
   errorVal[2] = errorAvg;
   errorVal[3] = errorPerOpAvg;
}

//~ copy(centers.begin(), centers.end(), ostream_iterator<Point3>(cout, "\n"));
